Overview

I know I billed this as “Ditching ArcGIS” in the poll, but I 1) am not an ArcGIS expert and 2) there may still be things that you must do in ArcGIS to get your work done. If that’s the case, there is the R - ArcGIS Bridge that makes it possible for these two analysis tools to play with one another, BUT I won’t discuss that here…

Configure chunk options

# Set global knitr chunk options
if (.Platform$OS.type == "unix") {
  # Do not specify Cairo device for MacOS
  knitr::opts_chunk$set(warning = F, message = F)
} else {
  knitr::opts_chunk$set(warning = F, message = F, dev.args = list(type = "cairo"))
}

Load libraries

You may need to install some packages (and their dependencies) that you don’t already have. The chunk below will do that for you automatically using the load_pkgs function.

# List packages required to run the script -------------------------------------
pkgs <- c("tidyverse","ggmap","sf","cowplot","here",
          "knitr","ggrepel","rnaturalearth","rgeos",
          "FedData","raster","mapview")
# Install and load all CRAN packages provided from a character vector
load_pkgs = function(pkgs) {
  new_pkgs = pkgs[!(pkgs %in% installed.packages()[ ,'Package'])]
  if (length(new_pkgs) > 0) install.packages(new_pkgs,repos = "http://cran.cnr.berkeley.edu/")
  invisible(lapply(pkgs,function(x)
    suppressPackageStartupMessages(library(x,character.only = T))))
}
# Load packages
load_pkgs(pkgs)

dir.create(here("Figs"))
dir.create(here("Output"))

A VERY basic map

# Load some ship track data
cps <- read_csv(here("Data/cps_nav.csv"))

# A basic 
plot(cps$long, cps$lat)

# A very basic ggplot2 scatter plot of the same data
ggplot(cps, aes(long, lat)) + geom_point()

A slightly less basic map

Let’s set a few parameters to make the map more spatially accurate. But it’s a long way from a useful “map”…

# Add some more aesthetics and coordinate controls, maybe some labels

# Read locations
locations <- read_csv(here("Data/locations.csv")) %>% 
  filter(group %in% c("city", "landmark") )

# Refine the "map"
ggplot(cps, aes(long, lat)) + 
  geom_point() +
  # geom_path(aes(group = transect, colour = factor(transect))) +
  # geom_point(aes(group = transect, colour = factor(transect), size = cps.nasc)) +
  # geom_point(data = locations, aes(long, lat)) +
  # geom_text(data = locations, aes(long, lat, label = name)) +
  coord_quickmap() +
  # coord_map(projection = 'azequalarea') +
  # coord_map(projection = 'azequalarea', xlim = c(-128,-123), ylim = c(45,50)) +
  theme_bw()

Canned mapping options

Before we go farther, let’s set some map boundaries based on either our data (cps$lat and cps.long) or manually (mb.lat and mb.long).

# Define lat and long bounds for west coast map
wc.lat  <- range(cps$lat)  #c(32, 52)
wc.long <- range(cps$long) #c(-130, -116)

# Define lat and long bounds for Monterey Bay
mb.lat  <- c(36.5, 37)
mb.long <- c(-122.2, -121.7)

# Set west coast boundaries for stamen maps
wc.bounds.stamen <- c(left = min(wc.long), bottom = min(wc.lat), 
               right = max(wc.long), top = max(wc.lat))
# Download stamen map of west coast
wc.map.stamen.toner <- get_stamenmap(wc.bounds.stamen, zoom = 6, maptype = "toner-lite") %>% 
  ggmap() + xlab("Longitude") + ylab("Latitude") + theme_bw()

# Add layers to map
wc.stamen1 <- wc.map.stamen.toner + 
  geom_point(data = locations, aes(long, lat)) +
  geom_text(data = locations, aes(long, lat, label = name)) + 
  geom_path(data = cps, aes(long, lat, group = transect, colour = factor(transect)),
            show.legend = F) +
  ggtitle("Basic options")

wc.stamen2 <- wc.map.stamen.toner + 
  geom_point(data = locations, aes(long, lat)) +
  geom_text(data = locations, aes(long, lat, label = name), 
            colour = "red", size = 2, hjust = 0, nudge_x = 0.5, angle = 45) + 
  geom_path(data = cps, aes(long, lat, group = transect, colour = factor(transect)),
            show.legend = F) +
  ggtitle("Formatted labels")

wc.stamen3 <- wc.map.stamen.toner + 
  geom_point(data = locations, aes(long, lat)) +
  geom_text_repel(data = locations, aes(long, lat, label = name), 
            size = 2, segment.colour = "black", segment.alpha = 0.5) +
  geom_path(data = cps, aes(long, lat, group = transect, colour = factor(transect)),
            show.legend = F) +
  ggtitle("Formatted and repelled lables")

# Combine maps
wc.grid.stamen <- plot_grid(wc.stamen1, wc.stamen2, wc.stamen3, nrow = 1)

# Save map
ggsave(wc.grid.stamen, filename = here("Figs/wc_map_stamen_toner.png"), width = 12, height = 7)
# Print map
include_graphics(here("Figs/wc_map_stamen_toner.png"))

# Set west coast boundaries for stamen maps
wc.bounds.google <- c(mean(wc.long),mean(wc.lat))
# Download stamen map of west coast
wc.map.google <- get_map(location = wc.bounds.google, zoom = 4, source = "google", maptype = "satellite") %>% 
  ggmap() + xlab("Longitude") + ylab("Latitude") + theme_bw() +
  xlim(wc.long) + ylim(wc.lat)

wc.google1 <- wc.map.google +
 geom_point(data = locations, aes(long, lat), colour = "white") +
  geom_text(data = locations, aes(long, lat, label = name), 
            colour = "white", size = 2, hjust = 0, nudge_x = 0.5, angle = 45) +
  ggtitle("Minor formatting")

wc.google2 <- wc.map.google +
 geom_point(data = locations, aes(long, lat), colour = "white") +
  geom_text_repel(data = locations, aes(long, lat, label = name), 
            colour = "white", size = 2) +
  ggtitle("Repelled labels")

# Arrange plots
wc.grid.google <- plot_grid(wc.google1, wc.google2, nrow = 1)

# Save map
ggsave(wc.grid.google, filename = here("Figs/wc_map_google.png"), width = 8, height = 7)
# Print map
include_graphics(here("Figs/wc_map_google.png"))

# Set west coast boundaries for stamen maps
mb.bounds.google <- c(mean(mb.long),mean(mb.lat))
# Download stamen map of west coast
mb.map.google <- get_map(location = mb.bounds.google, zoom = 10, 
                         source = "google", maptype = "satellite") %>% 
  ggmap() + xlab("Longitude") + ylab("Latitude") + theme_bw() +
  xlim(mb.long) + ylim(mb.lat)

mb.google1 <- mb.map.google +
 geom_point(data = filter(locations, name != "Monterey Bay"), aes(long, lat), colour = "white") +
  geom_text(data = filter(locations, name != "Monterey Bay"), aes(long, lat, label = name), 
            colour = "white", size = 4, vjust = 0, nudge_y = 0.01)
# Save map
ggsave(filename = here("Figs/mb_map_google.png"), width = 7, height = 7)
# Print map
include_graphics(here("Figs/mb_map_google.png"))

rnaturalearth

rnaturalearth for shoreline and country files. Below is a comparison of the medium- and large-scale country datasets.

# Import NaturalEarth coastlines and countries
coast_sf <- ne_coastline(scale = "medium", returnclass = "sf")
coast_sf_lg <- ne_coastline(scale = "large", returnclass = "sf")
countries_sf <- ne_countries(scale = "medium", returnclass = "sf")
countries_sf_lg <- ne_countries(scale = "large", returnclass = "sf")

# View unique regions
# sort(unique(countries_sf$region_wb))
# sort(unique(countries_sf$subregion))

# Select only certain regions to speed plotting
na_sf <- filter(countries_sf, subregion %in% c("Northern America","Central America")) 
na_sf_lg <- filter(countries_sf_lg, subregion %in% c("Northern America","Central America")) 

# Create West Coast map (medium scale)
wc.ne <- ggplot() + geom_sf(data = na_sf) +
  geom_point(data = locations, aes(long, lat)) +
  geom_text_repel(data = locations, aes(long, lat, label = name), 
                  size = 2, segment.colour = "black", segment.alpha = 0.5) +
  scale_x_continuous(name = "Longitude", limits = wc.long) + 
  scale_y_continuous(name = "Latitude", limits = wc.lat) +
  theme_bw() + coord_sf() +
  ggtitle("Medium scale")

# Create West Coast map (large scale)
wc.ne.lg <- ggplot() + geom_sf(data = na_sf_lg) +
  geom_point(data = locations, aes(long, lat)) +
  geom_text_repel(data = locations, aes(long, lat, label = name), 
                  size = 2, segment.colour = "black", segment.alpha = 0.5) +
  scale_x_continuous(name = "Longitude", limits = wc.long) + 
  scale_y_continuous(name = "Latitude", limits = wc.lat) +
  theme_bw() + coord_sf() +
  ggtitle("Large scale")

# Create map grid
wc.grid.ne <- plot_grid(wc.ne, wc.ne.lg, nrow = 1)

# Save map
ggsave(wc.grid.ne, filename = here("Figs/wc_map_natEarth.png"), width = 8, height = 7)
# Print map
include_graphics(here("Figs/wc_map_natEarth.png"))

Leaflet and mapview

mapview(ca_waters, alpha.regions = 0.1) +
  mapview(ca_mpas,zcol = "Type") +
  mapview(nasc.sf, cex = "bin.level", color = nascPal, alpha = 0.5)

GIS tasks using the sf package

Data packages

tigris

Access to the TIGER/U.S. Census line data files

# Get a map of Moss Landing
ML.bbox <- get_map("Moss Landing, CA", zoom = 14) 
# Plot the map
ggmap(ML.bbox) 
# Get the bounding box from its attributes
bb <- attr(ML.bbox, "bb")
# Create an extent polygon from the Monterey Bay bounding box
ML.extent <- polygon_from_extent(raster::extent(bb$ll.lon, bb$ur.lon,
                                            bb$ll.lat, bb$ur.lat), 
                             proj4string="+proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs")
# Map the bounding box to check the results
ggmap(ML.bbox) + geom_polygon(aes(x=long, y=lat), size=3, color='purple',
                             data=ML.extent, fill=NA)
# Download National Elevation Data
ned_ML <- get_ned(template = ML.extent, label="ned_ML", 
                    res="1", force.redo = F, raw.dir = here("Data"))
# Convert raster to data frame for plotting in ggplot
ned_df <- rasterToPoints(ned_ML) %>% 
  data.frame() %>% 
  rename(lat = y,
         long = x,
         elev = grdn37w122_1)

# Map elevation
ggplot(ned_df, aes(long,lat,fill = elev)) + geom_raster() +
  geom_contour(aes(z = elev), alpha = 0.5, binwidth = 1) + 
  scale_fill_viridis_c(name = "Elevation (m)") +
  theme_bw()

Other potentially useful packages

More R-GIS resources

Simple Feaures for R: Blogs, presentations, vignettes for the sf package

A Tidy Approach to Spatial Data: Simple Features from our own Eric Anderson

GIS with R: an introduction by Francisco Rodriguez-Sanchez (Pakillo)

An Intro to GIS with R by Jess Sadler (covers sp, sf, and rnaturalearth, among other things).